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The block LSMR algorithm for solving 
linear systems with multiple 
right-hand sides 


F. Toutounian* and M. Mojarrab 


Abstract 


LSMR (Least Squares Minimal Residual) is an iterative method for the 
solution of the linear system of equations and least-squares problems. This 
paper presents a block version of the LSMR algorithm for solving linear sys- 
tems with multiple right-hand sides. The new algorithm is based on the block 
bidiagonalization and derived by minimizing the Frobenius norm of the resid- 
ual matrix of normal equations. In addition, the convergence of the proposed 
algorithm is discussed. In practice, it is also observed that the Frobenius 
norm of the residual matrix decreases monotonically. Finally, numerical ex- 
periments from real applications are employed to verify the effectiveness of 
the presented method. 


Keywords: LSMR method; Bidiagonalization; Block methods; Iterative 
methods; Multiple right-hand sides. 


1 Introduction 


This paper is concerned with the solution of linear system of the form 


AX =B, AER", BER", s<n. (1) 
If A is large and sparse or sometimes not readily available, then iterative 
solvers may become the only choice. These solvers are categorized to the 
following three classes: 
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The first class is the global methods. The term global is due to Saad [34] 
and has been further expanded by Jbilou et al. [21] with the global FOM and 
GMRES algorithms for matrix equations. These methods are based on the 
use of a global projection process onto a matrix Krylov subspace. References 
on this class include [2,7,8, 12,13, 13, 21-23, 25-27, 32, 33]. 


The second class is the seed methods. The main idea of this kind of 
methods is briefed below. We first select a single system as the seed system 
and generate the corresponding Krylov subspace. Then we project all the 
residuals of the other linear systems onto the same Krylov subspace to find 
new approximate solutions as initial approximations. See [3,5,7,18,20,30,35] 
for details. 


The last class is the block methods which are more suitable for dense 
systems with preconditioner. The first block solvers are the block conjugate 
gradient (BI-CG) algorithm and the block biconjugate gradient (BI-BCG) al- 
gorithm proposed in [28]. Variable BI-CG algorithms for symmetric positive 
definite problems are implemented on parallel computers [19,29]. If the ma- 
trix is symmetric, an adaptive block Lanczos algorithm and a block version of 
Minres method are devised in [17]. For nonsymmetric problems, the BI-BCG 
algorithm [6,28], the block generalized minimal residual (BI-GMRES) algo- 
rithm [1,1,4,7,9-11, 36,37], the block quasi minimum residual (BIL-QMR) al- 
gorithm [14], the block BiCGStab (BI-BICGSTAB) algorithm [31], the block 
Lanczos method [34] and the block least squares (BI-LSQR) algorithm [15] 
have been developed. 


In this paper, we present a block version of LSMR algorithm [4] for solving 
the problem (1). Our algorithm is based on the block bidiagonalization [9]. 
We construct a simple recurrence formula for generating the sequences of ap- 
proximations {X;,} such that the Frobenius norm of A? Ry, decreases mono- 
tonically, where Ry = B— AXg. 


Throughout this paper, we use the following notations. For two n x s 
matrices X and Y, we define the following inner product: (X,Y) = tr(X7Y), 
where tr(Z) denoted the trace of the square matrix Z. The associated norm 
is the Frobenius norm denoted by || - || 7. We will use the notation (-,-)2 for 
the usual inner product in R” and the associated norm denoted by || - ||2. 
Finally, 0, and J, will denote the zero and the identity matrices in R***. 


The remainder of this paper is organized as follows. In Section 2, we give 
a sketch of the LSMR method and its properties. In Section 3, we present 
the block version of the LSMR algorithm. In Section 4, the convergence of 
the presented algorithm is considered. In Section 5, some numerical experi- 
ments on test matrices from the University of Florida Sparse Matrix Collec- 
tion(Davis [7]) are presented to show the efficiency of the method. Finally, 
we make some concluding remarks in Section 6. 
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2 The LSMR algorithm 


In this section, we present a brief of the LSMR algorithm [4], which is an 
iterative method for solving real linear system of the form 


Ax = b, 


where A is a matrix of order n and x,b € R”. 


LSMR. algorithm uses an algorithm of Golub and Kahan [10], which is 
stated as procedure Bidiag 1 in [32] to reduce the augmented matrix [b A] to 
the upper-diagonal form [61e, By], where e; denotes the first column of the 
identity matrix. The procedure Bidiag 1 can be described as follows. 


Bidiag 1 (Starting vector b; reduction to lower bidiagonal form) 
Au =6b, onv= ATu, 


BipiUi4 = Av; — aru, a 
a (De (2) 
O441Vi41 = AX Ui41 — Bi41%i, 


The scalars a; > 0 and 6; > 0 are chosen so that ||u;||, = ||v;||. = 1. With 
the definitions 


Q1 
Bz ag 
Uy, = [ui, ua, .-. Url, Vi, = [vi, v2, .-- UK]; B,= Tee . 
Be Ok 
Bri 
Lett = [Be Ontiertil, Veri = [Ve veri], 
the recurrence relations (2) may be rewritten as 
Un+i(Bi1e1) = 6, 
AV; = Ux4iBr, 
AT Un41 = VeBe + On410n41eh41 = Vert Lpat- 
Br 
AT AV, = ATU p41 Be = Vii Lig Be = Vert Re | Big, 
OR+1ER 41 
BE By 
= \V; K : 
eet Naess | 


This is equivalent to what would be generated by the symmetric Lanczos pro- 
cess with matrix A’ A and starting vector A7b. As we observe the procedure 
Bidiag1 will be stop if Av; — aju; = 0 or AP uj41 — i410; = 0, for some i. 
In exact arithmetic, we have UL Ue = I and VeEVE = I, where I is the 
identity matrix. 
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Hence using procedure Bidiag 1 the LSMR method constructs an approxi- 
mation solution of the form x, = Viyz which solves the least-squares problem 
min,, ||A7r,|], where r, = b— Ax. The main steps of the LSMR algorithm 
can be summarized as follows. 


Algorithm 1 LSMR algorithm 

Set Byuy = b, a1U1 = ATu, Q1 = 4, C1 = a1 61, Po 1, Po 1, Co 1, 

SO 0, hy U1; ho 0, Xo 0, 

For k = 1,2,..., until convergence Do: 
Proitet1 = AVE — ARUR, 
O641Vk41 = AT Un+1 — BaoiVe 
Pk = (aX + cPanen an 
Ck = Ak/pr, 

Sk = Br+i/Pk, 
Prot = SkAR41; 
Ore+1 = CkhOk+1; 
Ox = Sk—-1Pk, 


= .. 2 
Pr = ((Ck-1Pk)” + p41) ) 
Ck = Ck—-1Pk/Prs 


NI 


Sk = 941/Pr: 
Ck = CeCe, _ 
Crt = —SkCr» 


he = he — Oxpr/(Pr—1Pp—1))he—1, 

Te = Te-1 + (Ce /(PrPp) he, 

hepa = Ve+1 — (On41/pr)he, 

If |¢,.41| is small enough then stop, 
End Do. 


More details about the LSMR algorithm can be found in [4]. 


3 The block LSMR method 


We first recall the block Bidiag 1 algorithm [9]. This algorithm is the basis 
for our block LSMR method. 


The block Bidiag 1 procedure constructs the sets of the n x s block vectors 
Vi, Vo,... and U,,U2,... such that V,7V; = 0,, UFU; = 0s, for i # j, and 
VV; =F; UZU; = I,; and they form the orthonormal basis of R"***. 


Block Bidiag 1 (Starting matrix B; reduction to block lower bidiagonal 
form) 
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UB, =B, VijAy =A, 


Visi Bigi = AV; — U; AP, i-1.2 k 
Vist Ai = ATU — VBA, od see 


(3) 
where U;,V; € R°*°; By, A; € R°**, and U,B,,V, Ai, Ui41 Big1, Vier Aisi 


are thin QR decompositions of the matrices B, ATU;, AV; — U;A?, ATUj41 — 
VBA, respectively. With the definitions 


At 


By Ar 
Un = (Ui, U2, Aa A 


ee Ux], Ve = Mi, Vo, one 


- Vil, Th = oe 
B, AP 
Best 
the recurrence relations (3) may be rewritten as: 
Up Ei Bi = B, 
AV, = Opsilr, 


ATU p41 = VieTe t+ Vert Antibes, 


where EF; is the (A + 1)s x s matrix which is zero except for the rows i to 
i +s, which are the s x s identity matrix. We have also VV = I, and 
fo, 
Uys U e411 = Le41)s, where I is the | x | identity matrix. We define 

Test = [Te Evi Api | , 
then 


ae we esi 
ATU p41 = Vigil, 


” _ ao. oop _ ie 
ATAV, = ATO gaiTh = Vier Lggi tk = View kp ke 
Ani Eg 4 
TT, 


ay k 
= View eee | 2) 


At iteration k we seek an approximate solution X;, of the form 


Xp = Vez, 
where Y; is an ks X s matrix. Let By, = A,By for all k. Since 
A? R, = A’ B— AT AX, 


= V, A,B, — AT AV: Yp, 
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we have 
ARE VB KV Te Th 219 
i a Anyi Epis Tk : 
— = yaad 
= Vea (FB -— Buk Yz), (6) 
k 


where FE; is the ks x s matrix, which is zero except for kth s rows, which are 
the s x s identity matrix. 
In the block LSMR algorithm, we would like to choose Y; € R**** which 


minimizes the Frobenius norm of A? R,. From (6), A’ Ry can be written as 
E,By — TET. Y x 
ee 
— Beri ky Ye 


ATR, = Vivi 5 (7) 


where EF is the matrix obtained from E, by deleting its last block row. But 
since the columns of the matrix Vz4, are orthonormal, it follows that: 


phi pot, 2 
By Bye ry, 


Moet = fo go 
AT Rell = | = ||E,Bi—Te Tr Yelle + || Brot Es Yellz- 


Bae ee 
— Bri ky Ve a 
(8) 
We now define the linear operators yz and wx as follows. 
For Y € R*s*s 
xu(Y) = Ty TRY, 
and 
Lan oie 
We(Y) = Bri k;, Y. 
Then the relation (8) can be expressed as 
AT Rell = Ixe(We) — Ea Bille + live (Ve)IlB- (9) 


Therefore, Y, minimizes the Frobenius norm of the quantity A’R, if and 
only if it satisfies the following linear matrix equation 


XE (xe(Ye) — £1Bi) + UF (We (Ye) = Os, (10) 


where the linear operators XE and wr are the transpose of the operators x, 
and wx, respectively. Therefore, (10) is also written as the following 


band T T 
(TET)? (TE TeV, — F1B1) + (BraiE; )? (Brtily, Ye) = Os. (11) 


Hence, Y; is given by ne 
aH Te by 


where 
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“~ $e eee 4 es sj ~— 
Te = (Te Th)? + BeBe Beil: F, = T7T,E, Bi. (12) 


We define the matrix T;, as follows: 


As By 
” TIT, Bz Az 
— a aes —— . ‘s mers 5 
Byyiky early By, 
Br Ax 
Bri 


where A; = A;A? + BE Biss, for7 = 1,2,...,k. Therefore 

“A pa; eee ee oes a ses 

TSF Tie Be= (Arey (Babi) Os 4206), (13) 
and the approximate solution of the system (1) is given by 


Xp = VeTy Fr. 


Suppose that using the QR decomposition [11], we obtain a unitary matrix 


Q,;, such that 


Q1 By 03 
a2 Bs 4 
a i ee ial. Sa 
sXxks Or—2 Bp—1 Ok 
Ar-1 Pe 
Ak 


where R; is upper triangular as shown and @j;, 6;, 6; are the s x s matrices. 
So, 
ey ae “1 
Xp = VE(R, Re) Fr. 


By setting 
Py =Vik, = [Pi Po... Pa], 
and eee aa r 
F,=R, Fe= [ei v2 ---¥5] ; 
we have 


Py = (Ve — Pe-29% — Pr-1B,)&q > 
Xp = Xp-1 + Pryr. (15) 


From (15) the residual R;, is given by 
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Ry = Re-1 — AP ee, (16) 


where AP; can be computed from the previous AP;,’s and AV; by the simple 
update 7 7 
AP, = (AV, — AP 20% — APx-1B,)Oy - 


In addition, as [4], we show that the ||R;||~ can be estimated by a simple 
formula. By transforming 7; to block upper-bidiagonal form using a QR 


factorization: a. — Orit with One = P, we LP we have 
Ry = B- AX, 
= U, By — AV EY 5 


= U p41 (EB — TrY x) 


nan 


= O41 Qh 41 (Qep1 EB a a Yx). 


Since the columns of the matrices Q;41 and U;,4, are orthonormal, we have 


e R 
|Rulle = OvesE.Bs — [8] Yale a7) 
With definitions 
OBB SR Pe RES A a) 


the following Lemma shows that we can estimate ||R,||7 from just the last 
two blocks of Q,41£1B, and the last block of R,Yz. 


Lemma 1. In (17) and (18), 8; =% fori =1,2,...,k—1. 


aa: 


Proof. The proof is similar to that of Lemma 3.1 in [4] (see [28]). 


For the Frobenius norm of A’ R;, by using Theorem 1 (in section 4), we 
can also obtain the following simple formula: 


A? Rell = AP Re-illz — lleell, with | A? Rolle = ||Bille = llvollr. 


Now we can summarize the above descriptions as the following algorithm. 
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Algorithm 2 Algorithm (BI-LSMR ) 
Set Xo = Onxs, 
Set Gp = 05, b-1 Os, bo Is, Co Os, d_-1 =0,, do = Is, 
Set P_1 = Po =0nxs, 
Compute U,B, = B, V,A, = ATU; (QR decomposition of B and ATU), 
Set By, = A,B, 
Set y_1 =03, yo = —FPi, 
Set || A’ Rolle = ||vollr, 
For k = 1,2,..., until convergence Do: 
Wy, = AV; — U, Az, 
UrsiBrii = Wr (QR decomposition of W;,), 
Ar = Ay AP + Br a Beet, 
Sk = ATUR41 — VeBhiis 
Vert Agsi = Sp (QR decomposition of S;), 
Brat = Angi Bri, 
Br = dy-2 By, et 
On = Ch-1 BK + dh—-1Ak, 
By = Gp—1 Be + be-1 Ak, 
Bx = de-2By, ae 
Compute an unitary matrix Q(Gx, bp, Ck, dy) such that 


eal Belo 
Ce dk} [Bear] [0]? 

pe = —O," (Oy Pr-2 + Ba Pet); 

Py, = (Vig — Pr—29% — Pr—1By,)&y 

Xp = Xp-1 + Peeve, 

Ry = Re-1 — APeer, 

| AT Rell = AT Reale — ll vel? 

If || A? Ry || is small enough then stop, 
End Do. 


The BI-LSMR algorithm will be break down at step k, if @, is singular. 


This happens when the matrix Ls | is not full rank. So the BLLSMR 
k+1 


algorithm will not break down at step k, if By41 is nonsingular. We will not 
treat the problem of breakdown in this paper and we also assume that the 
matrices B;’s produced by the BI-LSMR algorithm are nonsingular. 


We mention that, we can use the BI-LSMR algorithm for computing a 
matrix solution X to the problem 


minimize||AX—Bllp, AE€R™*", BeR™S, s<min{m,n}, 


where m > n or m <n. In Section 5, we present the results of the BLLSMR 
algorithm for this kind of problems. 
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4 The convergence of the BI-LSMR algorithm 


In this section, we aim at studying the convergence behavior of the BI-LSMR 
method. We first give the following lemmas. 


Lemma 2. Let P;,1 = 1,2,...,k, be the n x s auxiliary matrices produced 
by the BI-LSMR algorithm and Ry be the residual matrix associated with the 
approximate solution X,, of the matriz equation(1). Then, we have 


(AT AP,)? AT Ry = 0s. 


Proof. Using Py, = VRS and equation(4), we have 


A? AP, = AT AP, Ex 


= ATAV,R, Ex 
iT 
= Tr | =-1 
=V eh S| Ree, 19 
k+1 Bae k k ( ) 


From (19), and (7), we have 


BB, — TERY, 


=T—-T T T 
(APAP,)T(AT Re) = Eg Ry (TET, (Ben Ex)"|ViwaWen [or 
—Brsrik; Yr 


—T—-T ~ T T 
=E,R, (TPT,(£.B, — TETY;,) — (Besiky )’ Basil, Ye) 
= 0s. (from (11)) 


= sap 
We note that Vx41 is orthonormal, thus Vz44Veq1 = L(r41)s- 


Lemma 3. Let P;,i = 1,2,...,k, be the n x s auxiliary matrices produced 
by the BI-LSMR algorithm. Then we have the following property 


P? A? AA? AP, = I,. 


Proof. Using (19), (12), (13) and (14), we have 
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(AT AP,)" (AT AP,) = (Vi i RB) (V ree Ro E,) 
a oc tel Bae. 4 a i+1 Bagh: ‘i ¥ 
ir 
—T—-T FEN aa Dede oN ei 
=EDR," |r?r, B7,,B,||5° or | RE: 
i [TPT Buy Bi BiyiB; i 


Sb Re TTR 


Pen eee, od ee ep) |" Res | a= ag 
= E; R; [Re Desa Q; Q; F = | RE; 


et jon [tie Oksx. | Fe E; 


sxks 


ay ae 


Theorem 1. Let X, be the approximate solution of (1), obtained from the 
BI-LSMR algorithm. Then 


|A7 Relle < ||A* R-all, 
where R, = B-— AXy. 
Proof. From(16), we have 
A? Ry_1 = AT Ry + AT APy yy. 

Using Lemma 2, since A? R; and A’ AP, are orthogonal, we have 

A? Ryall = ||A* Balle + ||AT AP evel. 
Thus 

| AT Rall = |AT Re—ille — ||AT APA Ell. 
Using Lemma 3, we have 


| AT Rell = AT Re_ill% — [IVell?, 
| AT Rielle < AT Reale. 


Theorem 1 is helpful in showing that if ||~x||7 is not very small in each 
iteration of the BI-LSMR algorithm, then the BI-LSMR algorithm will be 
stopped after a finite number of iterations. Otherwise, it is possible to occur 
stagnation. In this case, we can apply a reliable preconditioner for the block 
linear system of equations (1). 
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5 Numerical examples 


In this section, we consider the system AX = B, where AE R™*", Be 
R™**, X € R"**, and we present numerical results for several matrices taken 
from the University of Florida Sparse Matrix Collection (Davis [7]). These 
matrices with their properties are shown in Table 1. Our implementation is 
done on MATLAB version 07 on a PC machine with 4 GB RAM. Moreover, 
for the initial guess Xo = On. and B = rand(m, s), where the function rand 
creates an m X s random matrix with the coefficients uniformly distributed 
in [0,1]. The stopping criteria is set to || A? Rg||r/||Relle < 107'° x |All. 

Diagonal scaling was applied to the columns of [A, B] to give a scaled 
problem AX = B, in which the columns of [A, B] have unit 2-norm. By 
scaling, the number of iterations of BLLSMR for convergence reduced satis- 
factorily. 

In Table 2, we give the ratio t(s)/t(1), for s = 5,10, 20, and 30, where t(s) 
is the CPU time for BI-LSMR algorithm and ¢(1) is the CPU time obtained 
when applying LSMR for one right-hand side linear system. Note that the 
time obtained by LSMR for one right-hand side depends on which right-hand 
was used. So, t(1) is the average of the times needed for the s right-hand 
sides using LSMR. The results of Table 2 show that the BI-LSMR algorithm is 
effective and less expensive than the LSMR algorithm, because the indicator 
t(s)/t(1) is less than s. 

To show that the Frobenius norm of residual matrix decreases monoton- 
ically, we display the convergence history in Figure 1 for the systems corre- 
sponding to the matrices of Table 2 and BI-LSMR algorithm. In this figure, 
the vertical axis and horizontal axis are the logarithm in base 10 of the Frobe- 
nius norm of residual matrix and the number of iterations to convergence, 
respectively. We observe that for all matrices the Frobenius norm of residual 
matrix decreases monotonically. 

We display the convergence history of BLLSMR and BI-LSQR in Figure 
2 for the system corresponding to the matrix LPnetlib/lp_pilot. Figure 3 
(left and right) shows both solvers reducing ||A?R,||-/||Rxl|e and ||Rel|e 
monotonically and similarly. 
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SUL[OPOUL MOATOSOL [IQ GFZ ggg oul POIL FOI pueuroys/qH 

SULOoUISUS [ROTO EE OT O888T ou O0cE = =00GE loozeaps/teg 

suoryenbe [eluoelepip [eed — GZE O8E7 ou 006 006 o06epd/teg 

suoryenba Teluoloyrp [eyed =p S8S7T Ou. 1962 196¢ T96c9Pd/tTeg 

wofqoid surmurersoid reeuryT == G99 G99 ou LTE G0% = g0zoSs-d]/qreugT 

woTqoid surumurersoid reouryT = PG Gleyp = oul 098r = TPP = pt dy /quyoug'T 

wo[qoid Surumuresiso0id reoury CVO LZELOL ou 9961 978 soreurdt/qrjeug'T 

wo[qoid Surmuretsoid Ieoury ~—-9 6G very ou SEI 9g apTped[/qryeugT 

suoryenbe yetuerepip Teyqeq — 6ST PrLL  so& 006 006 0€-0-18/GH 

sutoysAs Joyndutod UT Serphys UOTyeTNUIIG —- TOT Icv ou SIT SIT GT 1018 /GH 
Soeur Teorey) 6FT V8IG ou 089 089 10895/G@H 

MUVNHONAP ddV oY} Ul posn xtyeur osreds wiopuey =~ TTS. POTES8T Ou OOOVT = O00rT ndde/uowg 
UsIsop MoM) orMoIpeTy OPE SP86r ou 0967 —dO96T ceppe/muey 

our diostq pl zuu wiAS suuIn[Oo smo Ayrodoig\Xtyey 
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Table 2: Effectiveness of BLLSMR algorithm measured t(s)/t(1) 


Matrix\s 5 10 20 30 
Hamm /add32 0.47 0.95 3.07 5.39 
Simon/appu 1.24 1.89 3.21 5.13 
HB/fs6801 0.27 0.38 0.97 1.19 
HB/grel115 0.99 0.51 3.41 8.57 
HB/gr-30-30 1.55 1.72 2.05 2.53 
LPnetlib/lpadlittle 0.37 0.42 1.63 12.54 
LPnetlib/lp_maros 2.92 3.75 6.79 12.36 
LPnetlib/Ip_pilot 2.40 4.95 15.90 22.92 
LPnetlib/Ip_sc205 0.70 1.30 2.11 4.70 
Bai/pde2961 0.33 0.52 0.98 1.14 
Bai/pde900 0.49 0.72 1.10 1.47 
Bai/rdb32001 0.30 0.39 0.38 0.76 
HB/sherman4 0.37 0.50 0.54 1.03 


=6= |p _adiittle 


‘= 0 = Ip_maros 


109 ollPylle 


Figure 1: Convergence history of the BIL-LSMR algorithm with s=20 


= = = BI-LSQR 
0 === BI-LSMR 


mm mB-LSQR 
——= BI-LSMR 


109, (IAT Ryl/llFyl) 
ke 
109 lll 


ee Camere ay Ses Pn 
Noo & k & Bw 4 © = b&w 


Figure 2: BL-LSMR and BI-LSQR solving a linear system AX = B with s = 20: problem 
LPnetlib/lp_pilot 
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6 Conclusion 


In this paper, we have presented a block version of LSMR algorithm for 
solving linear systems with multiple right-hand sides. We derived a sim- 
ple recurrence formula for generating the sequence of approximate solutions 
{X;,} such that the Frobenius norm of the quantity A’ R; decreases mono- 
tonically. In addition, we studied the convergence of the presented method. 
Besides, we showed that in absence of the break down condition, the pre- 
sented algorithm always converges. Numerical results have shown that the 
new algorithm obtains the results which are effective and less expensive than 
the LSMR algorithm applied to each right-hand side. 
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